## Loading required package: assertthat

## Warning: package 'assertthat' was built under R version 4.0.3

## 
## Attaching package: 'assertthat'

## The following object is masked from 'package:tibble':
## 
##     has_name

## Loading required package: glue

## Warning: package 'glue' was built under R version 4.0.3

## 
## Attaching package: 'glue'

## The following object is masked from 'package:dplyr':
## 
##     collapse

## Loading required package: rlang

## Warning: package 'rlang' was built under R version 4.0.3

## 
## Attaching package: 'rlang'

## The following object is masked from 'package:assertthat':
## 
##     has_name

## The following objects are masked from 'package:purrr':
## 
##     %@%, as_function, flatten, flatten_chr, flatten_dbl, flatten_int,
##     flatten_lgl, flatten_raw, invoke, list_along, modify, prepend,
##     splice

## `summarise()` has grouped output by 'ifbl_4by4', 'Jaar', 'TaxonIDParent'. You can override using the `.groups` argument.

## 
## -- Column specification --------------------------------------------------------
## cols(
##   id = col_double(),
##   code = col_character(),
##   naam_nederlands = col_character(),
##   naam_wetenschappelijk = col_character(),
##   taxn_vrs_key = col_character(),
##   taxongroep = col_character(),
##   usageKey = col_double(),
##   scientificName = col_character(),
##   rank = col_character(),
##   order = col_character(),
##   matchType = col_character(),
##   phylum = col_character(),
##   kingdom = col_character(),
##   genus = col_character(),
##   class = col_character(),
##   confidence = col_double(),
##   synonym = col_logical(),
##   status = col_character(),
##   family = col_character()
## )

## i Using ',' as decimal and '.' as grouping mark. Use `read_delim()` for more control.

## 
## -- Column specification --------------------------------------------------------
## cols(
##   substraat = col_character(),
##   substrate = col_character()
## )

1 Exploratory data analysis

Anthocerotophyta = hornworth = hauwmossen
Bryophyta = bryophytes = bladmossen
Marchantiophyta = liverworths = levermossen

phylum	n_obs	n_soorten	n_hokken
Anthocerotophyta	218	4	123
Bryophyta	95624	410	796
Marchantiophyta	17486	109	727

## `summarise()` has grouped output by 'phylum'. You can override using the `.groups` argument.

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

## Joining, by = "substraat"

phylum	substrate	number of species
Bryophyta	artificial stone	84
Bryophyta	bryophyte	3
Bryophyta	decaying vegetation	21
Bryophyta	decorticated wood	48
Bryophyta	epiphytic on living wood	73
Bryophyta	floating on water	1
Bryophyta	gravel or sand	108
Bryophyta	mineral soil	211
Bryophyta	peat	47
Bryophyta	rock, hard	95
Bryophyta	rock, soft	38
Bryophyta	soil on rock	68
Marchantiophyta	artificial stone	9
Marchantiophyta	bryophyte	17
Marchantiophyta	decaying vegetation	18
Marchantiophyta	decorticated wood	16
Marchantiophyta	epiphytic on living wood	15
Marchantiophyta	floating on water	3
Marchantiophyta	gravel or sand	34
Marchantiophyta	mineral soil	67
Marchantiophyta	peat	30
Marchantiophyta	rock, hard	14
Marchantiophyta	rock, soft	14
Marchantiophyta	soil on rock	20

phylum	ell_l	number of species
Bryophyta	1	9
Bryophyta	2	19
Bryophyta	3	37
Bryophyta	4	86
Bryophyta	5	125
Bryophyta	6	188
Bryophyta	7	189
Bryophyta	8	129
Bryophyta	9	15
Marchantiophyta	1	1
Marchantiophyta	3	20
Marchantiophyta	4	51
Marchantiophyta	5	51
Marchantiophyta	6	60
Marchantiophyta	7	56
Marchantiophyta	8	18

phylum	ell_r	number of species
Bryophyta	1	10
Bryophyta	2	63
Bryophyta	3	88
Bryophyta	4	67
Bryophyta	5	95
Bryophyta	6	141
Bryophyta	7	215
Bryophyta	8	101
Bryophyta	9	17
Marchantiophyta	1	28
Marchantiophyta	2	60
Marchantiophyta	3	20
Marchantiophyta	4	52
Marchantiophyta	5	30
Marchantiophyta	6	38
Marchantiophyta	7	21
Marchantiophyta	8	6
Marchantiophyta	9	2

phylum	ell_f	number of species
Bryophyta	1	12
Bryophyta	2	12
Bryophyta	3	61
Bryophyta	4	152
Bryophyta	5	187
Bryophyta	6	176
Bryophyta	7	51
Bryophyta	8	69
Bryophyta	9	54
Bryophyta	10	14
Bryophyta	11	4
Bryophyta	12	5
Marchantiophyta	4	15
Marchantiophyta	5	51
Marchantiophyta	6	52
Marchantiophyta	7	50
Marchantiophyta	8	47
Marchantiophyta	9	32
Marchantiophyta	10	9
Marchantiophyta	11	1

phylum	ell_n	number of species
Bryophyta	1	32
Bryophyta	2	195
Bryophyta	3	139
Bryophyta	4	196
Bryophyta	5	137
Bryophyta	6	71
Bryophyta	7	27
Marchantiophyta	1	43
Marchantiophyta	2	97
Marchantiophyta	3	45
Marchantiophyta	4	36
Marchantiophyta	5	23
Marchantiophyta	6	8
Marchantiophyta	7	5

phylum	ell_t	number of species
Bryophyta	2	52
Bryophyta	3	269
Bryophyta	4	142
Bryophyta	5	155
Bryophyta	6	108
Bryophyta	7	58
Bryophyta	8	13
Marchantiophyta	2	27
Marchantiophyta	3	113
Marchantiophyta	4	51
Marchantiophyta	5	43
Marchantiophyta	6	13
Marchantiophyta	7	2
Marchantiophyta	8	8

Only remove substrate levels with less than 5 species (for Ellenberg variables this is not necessary because they are included as continuous covariate).

## Registered S3 method overwritten by 'GGally':
##   method from   
##   +.gg   ggplot2

phylum	substrate	n
Bryophyta	artificial stone	84
Bryophyta	decaying vegetation	21
Bryophyta	decorticated wood	48
Bryophyta	epiphytic on living wood	73
Bryophyta	gravel or sand	108
Bryophyta	mineral soil	211
Bryophyta	peat	47
Bryophyta	rock, hard	95
Bryophyta	rock, soft	38
Bryophyta	soil on rock	68
Marchantiophyta	artificial stone	9
Marchantiophyta	bryophyte	17
Marchantiophyta	decaying vegetation	18
Marchantiophyta	decorticated wood	16
Marchantiophyta	epiphytic on living wood	15
Marchantiophyta	gravel or sand	34
Marchantiophyta	mineral soil	67
Marchantiophyta	peat	30
Marchantiophyta	rock, hard	14
Marchantiophyta	rock, soft	14
Marchantiophyta	soil on rock	20

## Warning: Ignoring unknown aesthetics: text

## `geom_smooth()` using formula 'y ~ x'

## Warning: Ignoring unknown aesthetics: text

## Warning in qlogis(prop_2000_2019): NaNs produced

## Warning in qlogis(prop_2000_2019): NaNs produced

## Warning in qlogis(prop_2000_2019): NaNs produced

## `geom_smooth()` using formula 'y ~ x'

## Warning: Removed 7 rows containing non-finite values (stat_smooth).

## # A tibble: 126 x 2
##    n_hok_1980_1999     n
##              <int> <int>
##  1               0   245
##  2               5    34
##  3               6    18
##  4               7    18
##  5               8    12
##  6               9    22
##  7              10    12
##  8              11    14
##  9              12    11
## 10              13     8
## # ... with 116 more rows

## # A tibble: 142 x 2
##    n_hok_2000_2019     n
##              <int> <int>
##  1               0    18
##  2               1    87
##  3               2    60
##  4               3    35
##  5               4    32
##  6               5    26
##  7               6    27
##  8               7    12
##  9               8    17
## 10               9    25
## # ... with 132 more rows

2 Methods

2.1 Trends

To assess the changes in distribution of bryophyte species form 1980 until 2018 we divided the survey in two periods and compared the period 1980-1999 with the period 2000-2018. To minimize the effect of badly prospected grid cells we only used grid cells (16km²) where at least 40 species were found in both periods (Figure 2.1).

16 km² grids of Flanders and the Brussels Capital Region with at least 40 species of Bryophytes recorded in both the period 1980-1999 and the period 2000-2019.

Figure 2.1: 16 km² grids of Flanders and the Brussels Capital Region with at least 40 species of Bryophytes recorded in both the period 1980-1999 and the period 2000-2019.

This resulted in 224 grid cells prospected in both periods. In these selected grid cells 35648 bryophyte observations (expressed as unique species per year and per 16 km² grid) were done in the first period and 45236 in the second period. To correct for potential changes in survey effort within those grids we used a relative change index as proposed by Telfer, Preston, and Rothery (2002), with slight modifications.

This method assumes that during the two surveys recorders attempted to observe as many species as possible in each grid cell. The size of the distribution area of individual species in a given period, i.e., the number of 16 km² grid cells where the species was observed, was expressed as a proportion of the total number of grid cells studied during that period. To summarise the changes in distribution area for all species, we modeled the number of occupied 16 km² squares in the second survey (out of the total number of 16 km² squares considered, i.e. proportional data) as a function of the logit-transformed proportion in the first survey (baseline) and covariates (phylum, substrate and Ellenberg values). Substrate effects were modeled for each phylum separately as random effects coming from a normal distribution with mean zero and an estimated standard deviation. Because the counts were overdispersed (variance larger than Binomial distribution), a beta-binomial distribution was used to account for overdispersion and a logit-link was used to link the proportions to the linear predictor with the following model equation:

\[ g(\mu_{i,s2}) = (\beta_0 + \beta_1\log\frac{\mu_{i,s1}}{1-\mu_{i,s1}} + \beta_L\text{ell}_L + \dots)*\text{phylum}_k + \mathbf{b}_{o,j} + \mathbf{b}_{1,j}\text{phylum}_k \] where:

\[ \begin{aligned} \mu_{i,s2} \equiv \mathbf{E}(Y_{i,s2}|\text{trials},\mathbf{b})&& \text{expected value}\\ g(\mu_{i,s2}) = \log\frac{\mu_{i,s2}}{1-\mu_{i,s2}} && \text{logit-link} \\ Y_{i,s2} \sim \text{BetaBinomial}(\mu_{i,s2}, \phi_{i,s2}) && \text{Beta-Binomial distribution}\\ \begin{bmatrix}\mathbf{b}_{o,j}\\ \mathbf{b}_{1,j}\end{bmatrix}\sim N(0,\Sigma) && \text{random intercept and slope}\\ \Sigma = \begin{bmatrix}\sigma^2_{0,j}&0\\0&\sigma^2_{1,j}\end{bmatrix} && \text{variance-covariance matrix} \end{aligned} \]

Species that occurred in less than five 16 km² grid cells in the first period were omitted form the analyses. We also used only species that were present in both periods.

A Bayesian hierarchical modeling procedure was used to calculate the posterior distribution for these parameters (Bürkner 2017) and we assumed non-informative priors.

We focused model selection on the Ellenberg terms only. All other effects were kept in the model. First, interaction terms for Ellenberg-values with phylum for which the 95% confidence interval contained 0 were removed. Next, main Ellenberg terms for which the 95% confidence interval contained 0 were removed. The simplest model with the highest predictive accuracy based on leave-one-out cross validation (Vehtari, Gelman, and Gabry 2017) was selected.

After rearranging the model equation (for one phylum) and exponentiating the linear predictor, ignoring all random effects, we obtain the following odds-ratio:

\[ \frac{Odds_{\mu_{i,s2}}}{Odds^{\beta_1}_{\mu_{i,s1}}} = \exp(\beta_0 + \beta_L\text{ell}_L + \dots) \] which allows us to predict the effect of Ellenberg values on the relative increase or decrease of the Bryophytes.

In a similar way, substrate random effects are calculated from the posterior linear predictor, conditional on Ellenberg values fixed at their median values. Gelman, Hill, and Yajima (2012) point out that no correction for multiple comparisons are necessary when substrate effects are calculated in this way.

All analyses were done using the R language for statistical computing (R Core Team 2020).

## Start sampling
## Start sampling

## Warning: Some Pareto k diagnostic values are slightly high. See help('pareto-k-diagnostic') for details.

Table 2.1: Results of model comparisons. First row contains the model with the highest predictive accuracy. The best model contained an effect for Ellenberg moisture and a phylum specific effect for Ellenberg temperature in addition to the effects dealing with sampling effort and substrate.
	elpd_diff	se_diff	elpd_loo	se_elpd_loo	p_loo	se_p_loo	looic	se_looic
alt_telfer_3	0.00	0.00	-3316.11	31.41	22.11	3.34	6632.21	62.82
alt_telfer_2	-0.88	1.27	-3316.99	31.10	23.61	3.27	6633.98	62.20
alt_telfer	-2.67	1.52	-3318.78	31.16	26.16	3.75	6637.55	62.32

plot(loo1$loos$alt_telfer, label_points = TRUE)

plot(loo1$loos$alt_telfer_3, label_points = TRUE)

alt_telfer_3$data[745,]

##     n_hok_2000_2019 n_hok_tot_2000_2019
## 745               8                 228
##     qlogis(n_hok_1980_1999/n_hok_tot_1980_1999)          phylum ell_f ell_t
## 745                                   -1.029619 Marchantiophyta     4     3
##     n_hok_1980_1999 n_hok_tot_1980_1999                substrate
## 745              60                 228 epiphytic on living wood

analyse_data_telfer_agg[745,]

## # A tibble: 1 x 14
##   phylum       taxon_id_parent n_hok_1980_1999 n_hok_2000_2019 n_hok_tot_1980_1~
##   <chr>                  <dbl>           <int>           <int>             <int>
## 1 Marchantiop~             760              60               8               228
## # ... with 9 more variables: n_hok_tot_2000_2019 <int>, prop_baseline <dbl>,
## #   prop_2000_2019 <dbl>, ell_l <dbl>, ell_n <dbl>, ell_f <dbl>, ell_r <dbl>,
## #   ell_t <dbl>, substrate <chr>

taxa[taxa$id == 760, ]

## # A tibble: 1 x 15
##      id code  naam_nederlands naam_wetenschap~ naam_duits naam_engels naam_frans
##   <int> <chr> <chr>           <chr>            <chr>      <chr>       <chr>     
## 1   760 ptil~ Boomfranjemos   Ptilidium pulch~ <NA>       <NA>        <NA>      
## # ... with 8 more variables: creation_date <dttm>, creation_user <chr>,
## #   update_date <dttm>, update_user <chr>, taxn_vrs_key <chr>,
## #   taxon_groep_id <int>, naam <chr>, omschrijving <chr>

#Boomfranjemos looks like an outlying observation
#it is the only epiphytic liverworth on living wood that declines strongly:
analyse_data_telfer_agg %>%
  filter(substrate == "epiphytic on living wood",
         phylum == "Marchantiophyta")

## # A tibble: 8 x 14
##   phylum       taxon_id_parent n_hok_1980_1999 n_hok_2000_2019 n_hok_tot_1980_1~
##   <chr>                  <dbl>           <int>           <int>             <int>
## 1 Marchantiop~             450             142             213               228
## 2 Marchantiop~             558             170             146               228
## 3 Marchantiop~             559             219             218               228
## 4 Marchantiop~             591              93             204               228
## 5 Marchantiop~             600               5              13               228
## 6 Marchantiop~             760              60               8               228
## 7 Marchantiop~             778              99             193               228
## 8 Marchantiop~             842              14               4               228
## # ... with 9 more variables: n_hok_tot_2000_2019 <int>, prop_baseline <dbl>,
## #   prop_2000_2019 <dbl>, ell_l <dbl>, ell_n <dbl>, ell_f <dbl>, ell_r <dbl>,
## #   ell_t <dbl>, substrate <chr>

3 Results

## Warning in checkMatrixPackageVersion(): Package version inconsistency detected.
## TMB was built with Matrix version 1.3.2
## Current Matrix version is 1.2.18
## Please re-install 'TMB' from source using install.packages('TMB', type = 'source') or ask CRAN for a binary version of 'TMB' matching CRAN's 'Matrix' package

## Registered S3 method overwritten by 'broom.mixed':
##   method      from 
##   tidy.gamlss broom

## Warning in tidy.brmsfit(alt_telfer_3): some parameter names contain underscores:
## term naming may be unreliable!

effect	component	group	term	estimate	std.error	conf.low	conf.high
fixed	cond	NA	(Intercept)	0.4904610	0.1315195	0.2316668	0.7477269
fixed	cond	NA	qlogisn_hok_1980_1999Dn_hok_tot_1980_1999	0.8964563	0.0173138	0.8631410	0.9301981
fixed	cond	NA	phylumMarchantiophyta	-0.7365657	0.1734234	-1.0744809	-0.3895969
fixed	cond	NA	ell_f	-0.0594441	0.0141049	-0.0868474	-0.0316259
fixed	cond	NA	ell_t	0.0265169	0.0186195	-0.0099043	0.0630969
fixed	cond	NA	qlogisn_hok_1980_1999Dn_hok_tot_1980_1999:phylumMarchantiophyta	0.0336040	0.0396186	-0.0438984	0.1119242
fixed	cond	NA	phylumMarchantiophyta:ell_t	0.0925428	0.0381588	0.0181477	0.1674919
ran_pars	cond	substrate	sd__(Intercept)	0.1613464	0.0518312	0.0851389	0.2810960
ran_pars	cond	substrate	sd__phylumMarchantiophyta	0.1092055	0.0826154	0.0042748	0.3041457

## Adding missing grouping variables: `phylum`

## Note: Using an external vector in selections is ambiguous.
## i Use `all_of(ellenberg)` instead of `ellenberg` to silence this message.
## i See <https://tidyselect.r-lib.org/reference/faq-external-vector.html>.
## This message is displayed once per session.

## Adding missing grouping variables: `phylum`
## Adding missing grouping variables: `phylum`
## Adding missing grouping variables: `phylum`

References

Bürkner, Paul-Christian. 2017. “Brms : An R Package for Bayesian Multilevel Models Using Stan.” Journal of Statistical Software 80 (1): 1–28. https://doi.org/10.18637/jss.v080.i01.

Gelman, Andrew, Jennifer Hill, and Masanao Yajima. 2012. “Why We (Usually) Don’t Have to Worry About Multiple Comparisons.” Journal of Research on Educational Effectiveness 5 (2): 189–211. https://doi.org/10.1080/19345747.2011.618213.

R Core Team. 2020. R: A Language and Environment for Statistical Computing. Manual. Vienna, Austria: R Foundation for Statistical Computing.

Telfer, Mark G, C. D. Preston, and Peter Rothery. 2002. “A General Method for Measuring Relative Change in Range Size from Biological Atlas Data.” Biological Conservation 107 (1): 99–109. https://doi.org/10.1016/S0006-3207(02)00050-2.

Vehtari, Aki, Andrew Gelman, and Jonah Gabry. 2017. “Practical Bayesian Model Evaluation Using Leave-One-Out Cross-Validation and WAIC.” Statistics and Computing 27 (5): 1413–32. https://doi.org/10.1007/s11222-016-9696-4.

Analysis of moss change

Wouter Van Landuyt, Hans Van Calster

2021-04-06

1 Exploratory data analysis

2 Methods

2.1 Trends

3 Results

References